{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 24\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rolling paper\n",
    "\n",
    "We'll start by loading the units we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "radian = UNITS.radian\n",
    "m = UNITS.meter\n",
    "s = UNITS.second"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And creating a `Params` object with the system parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = Params(Rmin = 0.02 * m,\n",
    "                Rmax = 0.055 * m,\n",
    "                L = 47 * m,\n",
    "                omega = 10 * radian / s,\n",
    "                t_end = 130 * s,\n",
    "                dt = 1*s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function estimates the parameter `k`, which is the increase in the radius of the roll for each radian of rotation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_k(params):\n",
    "    \"\"\"Estimates the parameter `k`.\n",
    "    \n",
    "    params: Params with Rmin, Rmax, and L\n",
    "    \n",
    "    returns: k in meters per radian\n",
    "    \"\"\"\n",
    "    Rmin, Rmax, L = params.Rmin, params.Rmax, params.L\n",
    "    \n",
    "    Ravg = (Rmax + Rmin) / 2\n",
    "    Cavg = 2 * pi * Ravg\n",
    "    revs = L / Cavg\n",
    "    rads = 2 * pi * revs\n",
    "    k = (Rmax - Rmin) / rads\n",
    "    return k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, `make_system` takes a `Params` object and returns a `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Make a system object.\n",
    "    \n",
    "    params: Params with Rmin, Rmax, and L\n",
    "    \n",
    "    returns: System with init, k, and ts\n",
    "    \"\"\"\n",
    "    init = State(theta = 0 * radian,\n",
    "                 y = 0 * m,\n",
    "                 r = params.Rmin)\n",
    "    \n",
    "    k = estimate_k(params)\n",
    "\n",
    "    return System(params, init=init, k=k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing `make_system`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can write a slope function based on the differential equations\n",
    "\n",
    "$\\omega = \\frac{d\\theta}{dt} = 10$\n",
    "\n",
    "$\\frac{dy}{dt} = r \\frac{d\\theta}{dt}$\n",
    "\n",
    "$\\frac{dr}{dt} = k \\frac{d\\theta}{dt}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes the derivatives of the state variables.\n",
    "    \n",
    "    state: State object with theta, y, r\n",
    "    t: time\n",
    "    system: System object with r, k\n",
    "    \n",
    "    returns: sequence of derivatives\n",
    "    \"\"\"\n",
    "    theta, y, r = state\n",
    "    k, omega = system.k, system.omega\n",
    "    \n",
    "    dydt = r * omega\n",
    "    drdt = k * omega\n",
    "    \n",
    "    return omega, dydt, drdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing `slope_func`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use an event function to stop when `y=L`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(state, t, system):\n",
    "    \"\"\"Detects when we've rolled length `L`.\n",
    "    \n",
    "    state: State object with theta, y, r\n",
    "    t: time\n",
    "    system: System object with L\n",
    "    \n",
    "    returns: difference between `y` and `L`\n",
    "    \"\"\"\n",
    "    theta, y, r = state\n",
    "    \n",
    "    return y - system.L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And look at the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final value of `y` is 47 meters, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "unrolled = get_last_value(results.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final value of radius is `R_max`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "radius = get_last_value(results.r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total number of rotations is close to 200, which seems plausible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "radians = get_last_value(results.theta) \n",
    "rotations = magnitude(radians) / 2 / np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The elapsed time is about 2 minutes, which is also plausible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_final = get_last_label(results) * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting `theta`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_theta(results):\n",
    "    plot(results.theta, color='C0', label='theta')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Angle (rad)')\n",
    "    \n",
    "plot_theta(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting `y`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_y(results):\n",
    "    plot(results.y, color='C1', label='y')\n",
    "\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Length (m)')\n",
    "    \n",
    "plot_y(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting `r`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_r(results):\n",
    "    plot(results.r, color='C2', label='r')\n",
    "\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Radius (m)')\n",
    "    \n",
    "plot_r(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also see the relationship between `y` and `r`, which I derive analytically in the book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.r, results.y, color='C3')\n",
    "\n",
    "decorate(xlabel='Radius (m)',\n",
    "         ylabel='Length (m)',\n",
    "         legend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the figure from the book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_three(results):\n",
    "    subplot(3, 1, 1)\n",
    "    plot_theta(results)\n",
    "\n",
    "    subplot(3, 1, 2)\n",
    "    plot_y(results)\n",
    "\n",
    "    subplot(3, 1, 3)\n",
    "    plot_r(results)\n",
    "\n",
    "plot_three(results)\n",
    "savefig('figs/chap24-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Animation\n",
    "\n",
    "Here's a draw function that animates the results using `matplotlib` patches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.patches import Circle\n",
    "from matplotlib.patches import Arrow\n",
    "\n",
    "def draw_func(state, t):\n",
    "    # get radius in mm\n",
    "    theta, y, r = state\n",
    "    radius = r.magnitude * 1000\n",
    "    \n",
    "    # draw a circle with\n",
    "    circle = Circle([0, 0], radius, fill=True)\n",
    "    plt.gca().add_patch(circle)\n",
    "    \n",
    "    # draw an arrow to show rotation\n",
    "    dx, dy = pol2cart(theta, radius)\n",
    "    arrow = Arrow(0, 0, dx, dy)\n",
    "    plt.gca().add_patch(arrow)\n",
    "\n",
    "    # make the aspect ratio 1\n",
    "    plt.axis('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "animate(results, draw_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Run the simulation again with a smaller step size to smooth out the animation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "**Exercise:** Since we keep `omega` constant, the linear velocity of the paper increases with radius.  Use `gradient` to estimate the derivative of `results.y`.  What is the peak linear velocity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(dydt, label='dydt')\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Linear velocity (m/s)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose the peak velocity is the limit; that is, we can't move the paper any faster than that.\n",
    "\n",
    "Nevertheless, we might be able to speed up the process by keeping the linear velocity at the maximum all the time.\n",
    "\n",
    "Write a slope function that keeps the linear velocity, `dydt`, constant, and computes the angular velocity, `omega`, accordingly.\n",
    "\n",
    "Run the simulation and see how much faster we could finish rolling the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
